; docformat = 'rst'
;
;+
;
; :Purpose:
;   Calculate weighted average two time series
;
; :Inputs:
;   tA: (input, double) time array defining mid-point (unless /start_of_interval set) of data xA. 
;     Must be regularly spaced, unless an output time array (tOut) is defined.
;   xA: (input) data array
;   tB: (input, double) time array for data xB, can be irregularly spaced
;   xB: (input) data array
;
; :Keywords:
;   dxA: (optional, input, float/double) uncertainty in A
;   dxB: (optional, input, float/double) uncertainty in B
;   dxOut: (optional, output, float/double) uncertainty in output
;   
;   nA: (optional, integer) number of data points contributing to each element in array A
;   nB: (optional, integer) number of data points contributing to each element in array B
;   nOut: (optional, output) number of data points in each element of output array
;   
;   start_of_interval: (optional, boolean) set to true if tA or tOut defines 
;     the start of the time interval for xA
;   monthly: (optional, boolean) set to true if tA is a regular monthly array
;   tOut: (optional, input) array of centre points of output array
;   replace: (optional, input) replace overlapping system A measurements with system B, 
;     instead of calculating an average of both
;   dxOut: (optional, output) array containing output uncertainties
;   
; :Outputs:
;   Regular or monthly array of combined data
;
; :Example::
; 
;   
;
; :History:
; 	Written by: Matt Rigby, University of Bristol, Jul 5, 2013
;
;-
function mr_combine_measurements, tA, xA, tB, xB, tOut=tOut, $
  dxA=dxA, dxB=dxB, $
  nA=nA, nB=nB, $
  start_of_interval=start_of_interval, monthly=monthly,  replace=replace
  
  
  ;Check if there are two datasets
  if n_elements(tB) gt 0 and n_elements(xB) gt 0 then begin
    if total(finite(xB)) gt 0 then begin
      two_datasets=1
    endif else begin
      two_datasets=0
    endelse
  endif else begin
    two_datasets=0
  endelse


  ;Define defaults
  if ~keyword_set(nA) then begin
    nA=finite(xA)
  endif
  if ~keyword_set(nB) and two_datasets then begin
    nB=finite(xB)
  endif
  if ~keyword_set(nB) and ~two_datasets then begin
    nB=0
  endif


  ;Define output time
  if ~keyword_set(tOut) then begin
    tOut=tA
  endif

  ;If no data, return
  if total(nA) + total(nB) eq 0. then begin
    return, {t:tOut, x:replicate(!values.f_nan, n_elements(tOut)), $
      dx:replicate(!values.f_nan, n_elements(tOut)), n:replicate(0, n_elements(tOut))}
  endif

    
  ;Define weights
  if keyword_set(dxA) or keyword_set(dxB) then begin
    if ~(keyword_set(dxA) and keyword_set(dxB)) and two_datasets then begin
      message, "If you set the uncertainty in one data set, you must set it in both"
    endif
    wA=1./dxA^2
    if two_datasets then wB=1./dxB^2
  endif else begin
    wA=replicate(1., n_elements(tA))
    if two_datasets then wB=replicate(1., n_elements(tB))
  endelse

  ;Calculate decimal dates
  tADec=mr_timedec(tA)
  if two_datasets then tBDec=mr_timedec(tB)


  ;Calculate output time series
  case 1 of
    keyword_set(monthly): begin
      caldat, tOut[0], mS, dS, yS
      caldat, tOut[-1], mE, dE, yE
      tOut=timeGen(units='MONTH', start=julday(mS, 15, yS, 0), final=julday(mE, 15, yE, 0))
      caldat, tOut, mOut, dOut, yOut
      caldat, tA, mA, dA, yA
      if two_datasets then caldat, tB, mB, dB, yB
    end
    keyword_set(start_of_interval): begin
      tOutDec=mr_timedec(tOut)
      dTDec=(tOutDec[1] - tOutDec[0])
    end
    else: begin
      tOutDec=mr_timedec(tOut)
      dTDec=(tOutDec[1] - tOutDec[0])
      tOutDec=tOutDec - dTDec/2.
    end
  endcase


  ;Define output arrays
  timeSize=n_elements(tOut)
  xOut=replicate(!values.f_nan, timeSize)
  dxOut=replicate(!values.f_nan, timeSize)
  nOut=replicate(0, timeSize)

  if keyword_set(monthly) then begin

    ;Calculate monthly means
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
    for ti=0, timeSize-1 do begin
      
      whA=where(finite(xA) and mA eq mOut[ti] and yA eq yOut[ti], countA)
      if two_datasets then begin
        whB=where(finite(xB) and mB eq mOut[ti] and yB eq yOut[ti], countB)
      endif else begin
        countB=0
      endelse

      ;Set countA to zero if we want to replace A with B
      if countB gt 0 and keyword_set(replace) then begin
        countA=0
      endif

      ;Calculate average
      if countA + countB gt 0 then begin
        
        av_array=fltarr(countA + countB)
        weight_array=fltarr(countA + countB)
        n_array=intarr(countA + countB)
        
        if countA gt 0 then begin
          av_array[0:countA-1]=xA[whA]
          weight_array[0:countA-1]=wA[whA]
          n_array[0:countA-1]=nA[whA]
        endif
        if countB gt 0 then begin
          av_array[countA:*]=xB[whB]
          weight_array[countA:*]=wB[whB]
          n_array[countA:*]=nB[whB]
        endif

        xOut[ti]=total(av_Array*weight_array)/total(weight_array)
        dxOut[ti]=sqrt(1./total(weight_array))
        nOut[ti]=total(n_array)
              
      endif
    endfor

  endif else begin

    ;Calculate regularly-spaced averages
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

    ;Remove time elements which correspond to non-finite data
    whA=where(~finite(xA), countA)
    if countA gt 0 then begin
      tADec[whA]=!values.f_nan
    endif
    if two_datasets then begin
      whB=where(~finite(xB), countB)
      if countB gt 0 then begin
        tBDec[whB]=!values.f_nan
      endif
    endif

    ;Calculate array indices using historgram
    histA=histogram(tADec, min=min(tOUtDec), Max=max(tOutDec), nBins=timeSize, locations=locations, $
      reverse_indices=riA)
    if two_datasets then begin
      histB=histogram(tBDec, min=min(tOUtDec), Max=max(tOutDec), nBins=timeSize, locations=locations, $
        reverse_indices=riB)
    endif
  
    for ti=0, timeSize-1 do begin
      countA=riA[ti+1] - riA[ti]
      if two_datasets then begin
        countB=riB[ti+1] - riB[ti]
      endif else begin
        countB=0
      endelse

      ;Set countA to zero if we want to replace A with B
      if countB gt 0 and keyword_set(replace) then begin
        countA=0
      endif

      ;Calculate average      
      if countA + countB gt 0 then begin
        
        av_array=fltarr(countA + countB)
        weight_array=fltarr(countA + countB)
        n_array=intarr(countA + countB)

        if countA gt 0 then begin
          av_array[0:countA-1]=xA[riA[riA[ti] : riA[ti+1]-1]]
          weight_array[0:countA-1]=wA[riA[riA[ti] : riA[ti+1]-1]]
          n_array[0:countA-1]=nA[riA[riA[ti] : riA[ti+1]-1]]
        endif
        if countB gt 0 then begin
          av_array[countA:*]=xB[riB[riB[ti] : riB[ti+1]-1]]
          weight_array[countA:*]=wB[riB[riB[ti] : riB[ti+1]-1]]
          n_array[countA:*]=nB[riB[riB[ti] : riB[ti+1]-1]]
        endif
  
        xOut[ti]=total(av_Array*weight_array)/total(weight_array)
        dxOut[ti]=sqrt(1./total(weight_array))
        nOut[ti]=total(n_array)

      endif
      
    endfor
  
  endelse

  return, {t:tOut, x:xOut, dx:dxOut, n:nOut}

end